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ABSTRACT 
A  system  of  two  partial  differential  equations  represents  the 
transient  heat  transfer  behavior  of  compact  heat  exchanger  surfaces 
when  subjected  to  a  step  change  in  fluid  temperature.   A  solution  is 
presented  for  this  system  which  includes  the  effects  of  longitudinal 
thermal  heat  conduction.   Also  presented  are  the  solutions  for  the  two 
limiting  cases  of  zero  and  infinite  longitudinal  conduction.   The 
numerical  results  have  been  compared  to  those  of  C.  P.  Howard  indicating 
a  significant  decrease  in  computational  time  and  an  increase  in  accuracy 
of  results.   The  revised  curves  of  maximum  slope  of  fluid  temperature 
versus  NTU  should  be  of  practical  value  in  the  evaluation  of  heat- 
transfer  data  obtained  by  transient  testing  of  compact  heat  exchanger 
surfaces.   An  unusual  combination  of  mathematical  techniques  is  presented 
for  the  solution  of  a  boundary  value  problem  involving  partial  differen- 
tial equations.   The  solution  combines  the  application  of  Laplace  trans- 
formation  with  a  numerical  technique  developed  by  H.  HurwitzsJr.,  and 
P.  F.  Sweifel,  and  adapted  by  L.  A.  Schmittroth  for  the  inversion  of 
Laplace  transforms.   This  technique  greatly  expands  the  number  of  cases 
to  which  Laplace  transforms  may  be  successfully  applied. 
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NOMENCLATURE 

2 
A        total  heat  transfer  area,  ft  . 

r     2 

A         matrix  cross  sectional  area,  ft  . 
s 

C         step  increase  in  fluid  temperature,  =  v  -v  , 
c         specific  heat  of  fluid,  BTU/lb  °F. 

c         specific  heat  of  solid,  BTU/lb  °F. 

2 
h         heat  transfer  coefficient,  BTU/hr  ft   °F. 

K         thermal  conductivity  of  matrix,  BTU/hr  ft  °F. 
s 

JL  total  length  of  fluid  flow  path,  ft. 

j\  dimensionless  conduction  parameter,  =  K  A  /w.c., 

' v  r  S  S   f  f 

NTU       dimensionless  heat  transfer  units,  =  hA/w  c  . 

t         dimensionless  time  variable,  =  hA0  /W  c  . 

*         s  s 

Q  time,  hr. 

u  solid  temperature,  °F. 

u  reference  solid  temperature,  °F. 

v  fluid  temperature,  °F. 

v^  fluid  temperature  at  X  =  0,  t  =  0  +;  °F. 

v  reference  fluid  temperature,  °F. 

wf  fluid  mass  flow  rate,  Ib/hr. 

W         mass  of  matrix,  lb. 

s  ' 

X         distance  along  the  fluid  flow  path,  ft. 
x         dimensionless  position  variable,  =  X  £  . 
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1.   Introduction. 

The  "single  blow"  problem  refers  to  the  study  of  the  transient  heat 
transfer  behavior  of  compact  heat  exchanger  surfaces  for  the  purpose  of 
determining  their  heat  transfer  characteristics.   Interest  in  compact 
heat  exchangers  has  grown  with  the  development  of  gas  turbine  engines. 
A  regenerative  cycle  is  desirable  to  increase  the  efficiency  of  the  gas 
turbine.   Proper  design  of  a  compact  heat  exchanger  requires  knowledge 
of  the  heat  transfer  characteristics  of  the  various  materials  in  numer- 
ous geometric  configurations.   Until  recently  it  has  been  necessary  to 
design  and  build  the  complete  heat  exchanger  in  order  to  determine  its 
heat  transfer  characteristics  experimentally.   The  objectives  of  the 
single  blow  studies  are  to  determine  the  heat  transfer  characteristics 
from  small  test  sections  of  the  various  matrices  and  to  obtain  consider- 
able reduction  in  experimental  costs. 

The  aim  of  the  experimental  testing  is  to  obtain  the  exit  fluid  time- 
temperature  history  subsequent  to  a  step  change  in  the  entrance  fluid 
temperature.   The  experimental  history  must  be  compared  to  a  theoretical 
time -temperature  history  to  determine  NTU,  the  dimensionless  heat  trans- 
fer parameter.   The  heat  transfer  coefficient,  h,  or  the  Colburn  J  factor 
can  then  be  calculated. 

In  1950  Locke  |jT]  outlined  five  experimental  techniques.   One  of  the 
simpler  of  these  has  become  known  as  the"maximum  slope"  method.   This 
technique  compares  the  maximum  slope  of  the  experimental  time -temperature 
history  to  the  theoretical  maximum  slopes  for  various  values  of  NTU  and 
of  /\    ,  the  longitudinal  conduction  parameter.   With  a  cursory  inspection, 
one  might  not  appreciate  the  efficiency  of  this  single  point  curve  matching 

Numbers  in  square  brackets  refer  to  the  bibliography. 


technique.   One  maximum  slope  represents  a  complete  theoretical  or  a 
complete  experimental  curve.   Therefore,  it  is  unnecessary  to  calculate 
and  plot  vast  quantities  of  theoretical  time-temperature  histories. 
This  method  eliminates  the  laborious  and  inaccurate  matching  of  theoreti- 
cal and  experimental  time -temperature  histories,  previously  required  to 
determine  NTU.   However,  the  maximum  slope  technique  has  two  major  limita- 
tions.  The  first  problem  is  the  difficulty  in  analytically  solving  the 
system  of  governing  differential  equations  when  longitudinal  conductivi- 
ties other  than  zero  or  infinity  are  included.   The  second  problem  is  the 
magnification  of  experimental  error  due  to  curve  matching.   This  problem 
is  discussed  in  detail  later. 

Considerable  analytical  effort  has  been  made  to  solve  various  aspects 
of  the  single  blow  problem  since  1927  when  Nusselt  first  studied  it. 
Schumann  |JL4_J  developed  a  solution  in  1929  for  zero  conductivity  in  the 
direction  of  flow.   His  results  are  for  a  porous  medium  such  as  gravel  but 
have  been  extended  to  include  matrices  consisting  of  metal  balls,  wire 
screens,  or  continuous  materials  for  which  the  effects  of  longitudinal 
thermal  conduction  can  be  ignored.   However,  as  indicated  by  Mondt  U-1J  and 
Creswick  3   „  usually  when  the  matrix  is  constructed  of  a  continuous 
material  in  the  flow  direction,  the  effects  of  longitudinal  conduction 
must  be  considered.   Creswick  outlined  a  finite  difference  technique  to 
include  this  effect  but  his  work  was  not  complete  enough  to  cover  the  area 
encountered  in  experimental  testing.   Mondt  (~ll"l  obtained  a  closed  solu- 
tion for  the  limiting  case  of  (K   —  OO   .   Figure  1  indicates  how  drasti- 
cally these  two  cases  differ  and  shows  the  need  for  intermediate  curves. 
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Figure  1 

In  1948,  F.  E.  Romie,  et.  al.jl.2J  developed  a  system  of  partial 
differential  equations  governing  the  transient  heat  transfer  behavior  for 
hollow  circular  cylinders.   Creswick  arrived  at  the  same  set  of  equations 
for  parallel  plates.   By  a  finite  difference  technique,  he  numerically 
approximated  the  solution  of  the  system.   However,  difficulty  with  con- 
vergence was  sometimes  encountered.   In  1963,  C.  P.  Howard  fsj  ,  guided 
by  physical  considerations,  developed  an  alternate  form  of  the  finite- 
difference  equations,  and  he  was  able  to  determine  the  convergence  criteria 
and  avoid  the  problem  encountered  by  Creswick.   Howard  presented  a  table 
of  results  and  a  set  of  curves  of  NTU  versus  Maximum  Slope  for  various 
values  of  f{   .   The  table  covers  a  range  quite  adequate  for  most  experi- 
mental testing.   As  the  values  of  NTU  increased  beyond  60,  it  became  very 
difficult  to  obtain  solutions  with  the  finite  difference  technique  due  to 
the  large  computation  times  involved. 

The  objective  of  the  solution  presented  in  this  paper  is  to  verify  or 
improve  the  results  obtained  by  C.  P.  Howard  and  to  avoid  large  computation 
times  as  NTU  increases.   This  solution  not  only  avoids  the  difficulty  at 


large  values  of  NTU,  but  becomes  faster  as  NTU  increases  due  to  the 
corresponding  increase  in  the  time  at  which  the  maximum  slope  occurs. 

Referring  to  the  error  magnification  due  to  curve  matching,  C.  P. 
Howard  Ts"]  showed  in  Figure  3-A  of  his  Appendix  3,  that  the  error  due 
to  the  experimental  determination  of  maximum  slope  is  magnified  by  a 
factor  of  two  or  greater  when  the  value  of  NTU  is  determined  by  the 
maximum  slope  technique.   The  error  magnification  factor  is  defined  as 
the  ratio  of  percent  error  in  NTU  to  percent  error  in  experimental 
Maximum  Slope.   The  essential  features  of  the  error  factor  as  a  function 
of  NTU  are  indicated  in  Figure  2. 
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Figure  2 
This  figure  shows  that  for  the  special  case,  /"I  =  0,  the  lowest 
error  magnification  is  between  2.5  and  2  for  NTU  greater  than  five,   At 
NTU  near  two  there  is  a  peak  for  which  an  accurate  determination  of  NTU 
by  maximum  slope  matching  is  impossible,   When  p{      is  non-zero, the  peak 


shifts  to  slightly  higher  values  of  NTU  and  the  useful  range  above  NTU 
of  four  decreases  as  n,      increases.   For  example,  for  ^  =  .06,  an 
error  factor  of  3,  at  NTU  =  5,  is  doubled  at  NTU  =  28.   Therefore,  with 
this  technique,  one  is  forced  to  accept  much  larger  errors  as  NTU  and 

ri      increase.   To  avoid  the  inherent  magnification  of  error,  a  modifica- 
tion can  be  made  to  the  time-temperature  curve  mai  ching  technique  with 
the  aid  of  modern  computers.   When  f\      is  finite,  non-zero  it  is  manda- 
tory to  use  a  computer  to  calculate  the  theoretical  values  of  fluid 

2 
temperature,  given  NTU,  and  t  .   It  would  be  elementary  to  include  a 

least  square  error  curve  matching  technique  in  the  temperature  program. 
Then,  given  three  or  more  experimental  values  of  fluid  temperature  and 
their  corresponding  values  of  t,  the  computer  would  compare  the  calcu- 
lated theoretical  temperatures  at  the  given  t's  for  various  NTU's  and 
determine  the  NTU  with  the  least  square  error.   The  anticipated  error 
in  any  particular  theoretical  value  of  temperature  is  +  1%  or  less. 


2 
"t"  denotes  the  dimensior.less  time  parameter, 


2.    Mathematical  Technique. 

The  approach  to  be  presented  solved  the  system  of  partial  differ- 
ential equations  for  finite,  non-zero  /L    by  eliminating  the  variable 
of  solid  temperature.   The  result  is  a  third  order  partial  differen- 
tial equation  of  fluid  temperature  as  a  function  of  time  and  x  (the 
ratio  of  distance  along  the  matrix  to  the  total  length).   This  equa- 
tion and  the  necessary  boundary  conditions  were  transformed  with  re- 
spect to  time  by  Laplace  transformations  to  obtain  a  third  order  total 
differential  equation  with  respect  to  x  with  parameter  s.   The  equation 
was  then  solved  for  the  transformed  fluid  temperature.   The  transformed 
boundary  conditions  were  applied  to  determine  the  three  arbitrary  co- 
efficients which  were  functions  of  the  parameter  s.   Because  the  trans- 
formed solution  was  very  complicated,  numerical  inversion  was  used  to 
compute  the  inverse  Laplace  transform. 

A  Gaussian  quadrature  method,  developed  by  H.  Hurwitz,  Jr.  and  P.  F. 
Sweifel  [oj  for  Fourier  Transform  Integrals,  was  adapted  by  L.  A.  Schmitt- 
roth  [_13]  for  the  numerical  inversion  of  Laplace  Transforms.   This  techni- 
que avoids  the  difficulty  encountered  with  alternating,  slowly  converg- 
ing functions  and  proved  to  be  essential  for  the  successful  inversion  of 
this  solution.   For  the  two  limiting  cases  f(   =  0  and  j\  =  OO  ,  a  direct 
inversion  of  the  Laplace  solution  was  available. 


3.    Solution. 

The  governing  differential  equations  for  the  transient  heat  transfer 

behavior  including  a  finite,  non-negative,  longitudinal  conductivity  are 
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as    follows 


dlX 


-(r\r-u)+A    -2j 


■z 


(i)         °}±     -     r\r-u)+^L    -d_U 

<9t        v        y    ntu  ax 

„,         Qar 

()       ~3X    =  NTU(a-Ar) 

The  applicable  boundary  conditions  are  the  following: 

a.  V(x,0)  =0  d.   3U  (o,t)  =  0 

9X 

b.  2^L(x,0)  =0        e.   3U  (l,t)  =  0 

9x  ax 

c.  v(0,t)  =  C 

The  following  is  the  approach  used  to  solve  this  set  of  equations 
for  the  fluid  temperature,  v,  and  the  slope  of  fluid  temperature,    Qfu 
respectively.   In  a  similar  manner  the  solid  temperature,  u,  and  other 
desired  values  can  be  obtained. 

Solving  equation  2  for  u  results  in 

(2a)   u  =    _X_^  Iff.    +.  or 

NTU  ax 

The  partial  derivative  of  u  with  respect  to  t  and  the  second  partial  of 
u  with  respect  to  x  are  determined  from  2a,  and  then  substituted  into 
equation  1  to  give 

3  2 

,,x       O  ■        3.  Qd*.t)    ,     MTj    dATfat)    _  NTU  _2flr(x,0 

c }  ax3  9X2-  A    ax 

_  NTU2_9or(xJt)  _   njjj.    9nr(x.t) 
R      dt  71  ax^t 


3 
The  equations  were  developed  from  a  heat  balance  by  Creswick  j"3~| 


Each  term  in  equation  3  is  transformed  with  respect  to  t  by  Laplace 
transformation  along  with  the  following  boundary  conditions: 
a.   v(x,0)  -  0 

b.  2LQ£  (x,o)  =  o 

The  necessary  transforms  are 

etc 


oo  oo        ,<x>  si  \  i 

<^  cat         J      4  3t  o     'o 

=  -  or  Cx,  0)  \  5  [nr  ex,  S)J 

Substitution  of  these  values   in  equation  3  after  applying   the  above  bound- 
ary conditions,    results   in 

^      gATCX.S),  NTU  32AT(x,s)      NTU^S-H^/VCx.S)      njTu2S  Wfos) 

(3a)-°~o)x3  2/2         ^  3x  R 

The  corresponding  auxiliary  equation  is 


(3b).  n3  +  h  ru  n?  -  iffi^  (s+d)  n  -  ^nL  S  =■  O 


The  general  solution  in  the  Laplace  s-plane  for  the  fluid  temperature  is 

where  rl,  r2,  r3   are  the  roots  of  equation  (3b). 

The  boundary  conditions  C ,  d,  and  e  are  transformed  and  then  used  to 

determine  the  coefficients  CI,  C2,  and  C3.   The  transform  of  C   is 

s 


The   transform  of  d  is  iZ-LLCOjS)  =  Q      ,    and  similarly   for  e, 

3U.  (IS)  -0 
dX 

From  equation  2a  the  following  equation  is  obtained; 


(5).      Su^Cx.t)        J^_   9V       ,      9nr(x,t) 
<3x  ~    ntu   ax2  ax 


Applying  boundary  condition  *  to  equation  4  gives 


(6).    nr<:o,s)   =  Cl(S)  +  Cz(s)-h  £3cs)  =  C_  . 

s 


Rather  than  apply  boundary  conditions  d  and  e  directly  it  is  convenient  to 
use  their  equivalent  form  in  terms  of  v.   To  apply  boundary  condition  d 
it  is  necessary  to  transform  equation  5  which  results  in 

i£  (X,S)_  A__  £j2T(X,S)      $>nr()(,S)  therefore; 

ax         "  ntu  sx*        +  ax 

sx  ntu  3X2-  ax  • 

When  this  boundary  condition  is  applied  the  following  equation  is  obtained: 

U;*  NTU 
The  application  of  boundary  condition  e>  is  similar,  resulting  in 

(8).  Un^en^AlCzen^n.^eni)^^^  W^e*  Wie*^  0, 

NTU 

(8a).   Define   R^  =  f  ArD_  4-Hml   where  n  =  1,  2,  3.   Then  the 
following  matrix  equation  can  be  set  up  using  equations  h*    7,  8. 
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To  solve  for  CI,  C2,  C3,  Cramer's  Rule  is  applied.   The  denominator  is 


a=  [R,R3(enL  eft2  )+  R.R^'-e*1)  +  R,R,.(e'*_efc  J] . 


The  three  arbitrary  coefficients  are  the  following: 
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and 


R.Rse"'  -R,R3e*3 
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Substituting  CI,  C2,  and  C3  into  equation  4  results  in 
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5  /\ 


Evaluating  equation  9  at  x  =  X     results  in  the  following: 
(10) 


d,s)  !1  dgdi^^ 


("2+^3)_  AM$][ 


Recall  that  the  Laplace  operation  S-p(S)  —  P(+0)      corresponds  to 

iL  F(t).     It  follows  directly  that   -^^  (i>t)  is  equal  to 

dt  ^t 

the  operation     S-/ir(l,S)         »  since  /Vfi^)—    O      b? 
boundary  condition   cl   .   Therefore, 

To  determine  the  exit  fluid  temperature  and  the  slope  with  respect 
to  time  of  the  exit  fluid  temperature,  it  is  necessary  to  find  the  in- 
verse Laplace  transforms  of  equations  9  and  10  respectively.   If  the 
values  of  R.  and  r,,  i  =  1,  2,  3;  were  not  functions  of  s  this  step  would 
be  elementary.   Unfortunately,  the  values  of  R  and  r  are  indeed  functions 
of  s  as  indicated  by  equations  8a  and  3b  respectively.   Therefore,  for 
this  particular  case  the  task  of  finding  an  analytic  inverse  transform 
is  nearly  hopeless. 

Rather  than  proceed  further  on  this  approach,  a  computer  program  to 
calculate  the  complex  roots  for  given  values  of  s,  NTU,  and   ^  ,  is 
used.   Consequently,  this  step  precludes  any  possibility  of  finding  an 

analytical  solution  if  one  exists.   Now,  the  value  of  v(l,s)  or       at 

dt 

(l,s)  can  be  determined  for  any  set  of  the  parameters,  C,  s,  NTU,  and   A   , 
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and  the  Laplace  inversion  integral  may  be  used  to  find  v(l,t).   Find- 
ing an  adequate  numerical  technique  to  evaluate  this  inversion  integral 
was  no  simple  matter. 

When  Simpson's  Rule  was  applied  to  evaluate  the  inversion  integral 
using  the  real  (VR)  and  imaginary  (VI)  parts  of  v(l,s),  the  results  were 
grossly  in  error,  the  behavior  being  as  described  in  reference  6 1  and  Il3j . 
The  numerical  technique  described  in  reference £l 3] was  found  to  be  greatly 
superior  to  Simpson's  Rule  for  this  particular  problem.   It  incorporates 
a  technique  to  accelerate  convergence  of  an  alternating  series,  and  uses 
Gauss-Chebyshev   formulas  for  integration.   This  technique  uses  either 
the  real  part  (VR)  or  the  imagirary  part  (VI)  and  results  in  two  differ- 
ent forms.   Since  there  are  no  known  values  for  this  solution  it  was 
necessary  to  use  both  of  these  forms  to  be  confident  that  the  correct 

solutions  were  obtained. 

4- 

According  to  the  theory  of  Laplace  inversion,  the  inversion  integral, 

1        .      WCXiS)  P       nS  ,    should  be   independent   of  G,    where 

2TTL    Jgwoo  wc       Ui 

s  =  G  +  iy  .   When  this  approach  was  used  with  Simpson's  Rule  none  of 
the  results  were  independent  of  G.   Recall  the  form  of  the  inversion 

integral     _S —  J  /V"(XiGi  <j)  ^^  (j H  '      Since  the  re8ion 

initially  investigated  used  t  =  7  and  G  of  .5,  1  and  2,  any  error  in  the 
numerical  technique  was  greatly  multiplied  giving  very  poor  results,  However; 
this  indicated  that  G  is  an  additional  parameter  that  must  be  determined 
for  each  new  program.   Rather  than  cause  difficulty,  this  parameter  proves 
to  be  essential.   For  the  same  G,  it  was  found  that  one  program  gave  con- 
sistently low  results  while  the  other  was  consistently  high.   Then  by 


4 
Churchill,  op.  cit.,  p.  176. 


12 


varying  G  of  both  programs  the  program  using  VR  can  be  made  to  agree 
with  the  results  of  the  program  using  VI.   This  is  the  approach  that 
was  used  to  make  the  various  programs  converge  to  a  solution.   The 
numerical  technique  outlined  in  |_13J  uses  Chebyshev' s  polynomial  to  fit 
the  functior   f(s).     It  was  found  that  this  polynomial  is  very  sensi- 
tive to  a  small  change  in  G.   For  example,  upon  changing  G  from  0.0  to 
0.05,  the  product  (slope t  NTU)  changed  from  1.0877  to  1.1333  an  increase 
of  4.19%.   Again  according  to  theory,  results  should  be  independent  of  G. 
Therefore,  two  new  programs  were  written  using  a  Gauss-Legendre  quadra- 
ture format  to  integrate  VR  or  VI.   These  programs  were  found  to  be 
essentially  independent  of  G  for  reasonable  increments  of  G.   For  example 
changing  G  from  0.0  to  0.06,  the  slope»NTU  changed  from  1.10096  to  1.09776 
a  decrease  of  0.29%.   Figure  3  indicates  a  comparison  of  these  two  pro- 
grams.  On  this  basis  alone,  it  is  felt  that  the  Legendre  polynomial  is 
more  accurate.   Unfortunately,  in  comparison  to  the  programs  using  the 
Chebyshev  polynomial  it  is  considerably  slower.   Because  of  this  fact, 
the  Legendre  programs  which  include  an  error  analysis  are  being  used  to 
determine  accurate  test  values  for  the  determination  of  the  best  G  for 
each  program.   Once  the  best  G  is  determineds  the  speedier  Chebyshev 
programs  can  be  used  for  long  data  runs. 

The  integration  technique  described  in  detail  in  references  [6 J  and 
IJ-3J  will  be  described  here  as  it  was  applied  to  this  problem.   The  inver- 
sion integral  is  the  following: 

<12)    Z  Xm)=m  ==£-/.     P«)e-Stds 

Then  writing  f(s)  as  the  sum  of  the  real  and  imaginary  parts  results  in 
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The  inversion  integral  can  then  be  written  as 
(12a) 


W=~  yJftte+L$+W(G*^ 


where  the  following  condition  applies:   f(s)   =  f(s). 
It  follows  directly  that 


(12b) 
(12c) 


and 


In  the  first  integral  of  equation  12a,  let  y  =-y\    and  dy  =  -dy* ,  then 

Applying  the  relations  12b  and  12c,  equation  13  reduces  to 

St  ,oo  -n 

W=  2?  J0  [[* Cg+l^-^ WG-Ky)]eLy  +[^re^)+i- W^^>]ety J^ 

On  rearranging  and  simplifying,  this  becomes 


&"t  ,oo 


(i4>£ft)-iL_^  [0CG+cy)cosyt-LH^+^j5inyfJdLj^ 


By  definition  of  the  inverse  integral,  f(-t)  =  0.   Thus 


(15) 


-Gi^°9- 


A  bar  over  a  quantity  indicates  the  complex  conjugate, 
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which  reduces  to 

Two  forms  of  f(t)  can  be  found  by  expressing  the  integrals  either  in 
terms  of  the  real  or  the  imaginary  parts.   They  are  respectively 


(16) 


Ut)-^~     J      <P(G+Ly)C0Syt  dtj,     and 


TT  yo 

(17)      kt)  =--^^°Vrs«V)  Sinqt  Jy . 

These  two  equations  are  the  basis  of  the  two  programs  mentioned  earlier. 
In  general  discussions  VR  will  represent  the  function  d)  cos  yt  and  VI 
the  function  fp  sin  yt.   The  Gaussian  quadrature  formula  using  either 
Chebyshev's  or  Legendre's  polynomial  is  used  to  evaluate  the  integral 
from  zero  to  infinity   Both  functions  VR  and  VI  are  of  the  general 
form  shown  in  Figure  4 


-VR 

Or 


Figure  4 
If  the  function  were  integrated  with  respect  to  y  and  the  summation 
stopped  at  a,  there  would  be  significant  error.   Two  operations  are 


16 


included  in  this  solution  to  reduce  the  error.   The  first  operation 
is  an  averaging  technique  which  accelerates  convergence.   The  second 
operation  involves  integrating  each  half  wave  and  then  summing  them. 
To  accomplish  this  we  introduce  the  following  change  of  variable  in 
Equation  17. 

Let     X—  —  V  -  4-  therefore    U  -^(X+'Z )  and    Ju-~dJC% 

Applying  these  to  equation  17  results  in 

■oo 


(17a)    f(t)=-2£^)i    ^[G+^(x+x)JsinTT(x+^)dx, 

co    rk.J^tfx-    i —  — i 

Z.     )'        ^  G  +  L^Cx+^sinTr(x4)J/. 


t 


Now  let      x'—    X-'n     ,    then 


n    /o^t-      CO 


/>a  =  o   /-  /^ 


,  £       r  -, 


Upon  supressing  the  prime,  we  have 


(is)  f(t)=  2e  J(-_1),n«J(2z//[G4.i|:(x+m4J]cosTrxdx 

When  equation  18  is  written  in  terms  of  partial  sums,  the  result  is 


(19) 


HO 


2  pjk*     oo 
— r—   Z  1^  It) 


/>i  =  o 


where 


and 
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(19a)  Im(t)=(-1)     I,  ^ 


G  -H  L  51  C-X  +  m  +-i-)l  c  OS  cttx;  dx 


Again  referring  to  j  131  ,  L.  A.  Schmittroth  applied  a  Gaussian  quadrature 

method  to  evaluate  I  (t).  I  ,'t)  can  b   approximated  as  follows: 

n       n 


j,!  COSTAL  J 


where 


2(2  N  4-1) 


N) 


and  _  uj  -  TT 


t 


(-^♦^j^eHH-rf+m-Hs) 


The  number  of  points  used  to  fit  the  function  is  2  N,   Obviously,  the 
larger  N  is,  the  greater  the  degree  of  accuracy  one  can  expect.   This 
is  the  core  of  the  iterative  process  and  an  adequate  N  is  desired  to  en- 
sure reasonable  computer  time.   This  point  will  be  discussed  in  detail 

N 
later.   The  N  coefficients  W   are  the  solutions  of  the  system. 


KJ 


<19d)    2  Z  wNcos 


ZN-2 


J  =  * 


(2J-l)TT 
2C2N+1) 
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where      h   =   1,    2,    . „ . ,    N. 

Equation  19  uses  Chebyshev's  polynomial  for  the  curve  fitting.   A 
similar  form  of  T  t)  can  be  developed  from  Equation  16,   As  mentioned 
earlier,  an  integration  technique  was  developed  using  a  Legendre  poly- 
nomial for  curve  fitting.   If  equation  16  is  used  to  illustrate  this 
approach,  thr  relation  *  =  tv_   is  required.   Then  equation  16  reduces  to 
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CO 


Ct  r  l 

(20)  f(i)  =  l2-      0(G+l*:x)CosTrxdx. 

"to  ^ 


Now,  the  Gaussian  quadrature  method  can  be  applied  as  above.   For  the 
Gauss-Legendre  formulas,  a  change  of  variable  was  applied  to  equation  17b 
to  set  up  the  infinite  sum  of  integrals  from  zero  to  one,  one  to  two  and 
so  on  rather  than  a  sum  of  integrals  evaluated  from  -1/2  to  1/2. 

The  Gauss-Chebyshev  and  the  Gauss-Legendre  formulas  applied  to 
the  real  or  imaginary  form  of  the  inversion  integral  produce  four  differ- 
ent computer  programs.   Any  one  of  these  should  give  the  correct  solution, 
but  one  form  may  prove  to  be  better  for  a  particular  application. 

For  zero  and  infinite  longitudnal  conduction,,  the  governing  differ- 
ential equations  reduce  to  much  simpler  forms.   The  solutions  for  both 
cases  are  given  in  Appendix  I  and  are  in  complete  agreement  with  the 
solutions  of  Schumman  and  Mondt  for  r\   =  0  and  *f\    -   OQ    respectively. 
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4.   Computer  Program. 

The  principal  objective  of  the  computer  program  is  to  solve  for  the 

6  n 

maximum  slope  of  fluid  temperature  at  X  -  X*  .   The  main  program  is  design- 
ed to  determine  the  slope  of  the  fluid  temperature  at  X  -H  which  is  the 
numerical  inversion  of  equation  11. 

The  subroutine  S00TS2  calculates  the  real  and  imaginary  coefficients 
of  the  equation  3b  given  the  parameters  NTU,  r\    ,  G,  and  Y.   These  co- 
efficients are  then  put  into  subroutine  T00TS2. 

Subroutine  T00TS2  calculates  the  real  and  imaginary  parts  of  the 
roots  r  ,  r  ,  r  of  equation  3b.   It  sends  its  results  back  to  S00TS2 
which  arranges  them  in  the  proper  order  and  subscripts  them   appropriately. 
S00TS2  then  provides  the  proper  roots  as  input  values  to  EVAL.   The  EVAL 
subroutine  provides  the  main  program  with  the  values  of  the  real  (VR) 
and  imaginary  part  (VI)  of    £__-  ( J-}  5  )     ,  given  the  parameters  G, 

at 

IT  X,  C,  and  the  roots  from  S00TS2. 
t 

The  main  program  will  call  for  N  of  these  values  (see  equation  19c) 

with  which  it  will  calculate  I  (t)  (see  equation  19).   The  programmer  is 
free  to  determine  the  number  of  I  (t)'s  to  sum,  where  ^ .  ~J_rr\(t)~   3/yn('t)t 
Fifteen  to  twenty  partial  sums  appear  to  be  sufficient  for  the  function  en- 
countered in  this  problem.   Each  partial  sum  calculated  is  stored  and 
serves  as  an  input  to  subroutine  AVER. 

The  purpose  of  subroutine  AVER  is  to  accelerate  convergence  of  the 
partial  sums  by  an  averaging  technique.  This  technique  is  discussed  in 
reference  I  6  J  .   The  subroutine  will  take  an  arbitrary  number  of  partial 
sums  and  calculate  any  specified  n   average,  and  return  it  to  the  main 
program.   It  also  calculates  the  difference  between  successive  n   averages 

6i.e.,  9QT. 

dt 
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returning  it  to  the  main  program.   The  main  program  then  requires  this 
difference  to  be  less  than  an  arbitrarily  specif iec  value  (EP).   If  the 
difference  is  larger  than  EP,  the  main  program  increases  the  number  of 
partial  sums  by  LD  and  calls  AVER  again  until  the  difference  is  less 
than  EP.   At  this  point,  the  value  of  the  accelerated  partial  sum  is 
accepted  and  the  slo^e  is  calculated.   However,  the  slope  calculated  for 
a  particular  value  of  dimensionless  time  is  not  the  value  desired.   It 
is  necessary  to  search  for  the  maximum  slope  by  varying  time. 

Initially,  a  search  using  an  incremental  step  of  time  was  used. 
This  was  eventually  abandoned  because  it  would  find  relative  maxima  only 
and  it  was  tediously  slow  if  a  poor  choice  were  made  for  a  starting  time. 
It  should  be  noted  that  theoretical  relative  maxima  have  been  found  for 
small  values  of  NTU  and  negative  slopes  have  also  been  noted.   J.  M. 
Bannon|_l|  is  concurrently  conducting  an  experimental  thesis  and  has  experi- 
mentally found  relative  maxima  and  negative  slopes  for  high  mass  flow  or 
small  values  of  NTU.   It  would  be  interesting  to  compare  a  theoretical 
time  temperature  curve  for  his  experimental  value  of   J\  with  his  ex- 
perimental time -temperature  curve. 

The  purpose  of  SLOMAX  is  to  find  the  maximum  slope.   To  avoid  the 
problem  of  choosing  a  poor  starting  time,  a  reasonable  range  of  t   is 
chosen.   A  rule  of  thumb  is  that  the  largest  maximum  slope  will  occur  at 
a  time  less  than  or  equal  to  NTU.   SLOMAX  uses  a  parabola  curve  fitting 
technique  with  three  points. 

The  first  and  last  of  the  three  points  are  program  input  values  of 
time  and  the  middle  point  is  the  arithmetic  mean.   The  slopes  are  calcu- 
lated for  the  first  two  values  of  time  and  then  tested  to  ensure  that  the 
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middle  slope  is  larger  than  the  first.   This  is  necessary  to  obtain  the 
maximum  rather  than  a  minimum  slope.   Then  the  third  value  of  slope  is 
calculated  and  tested  to  insure  that  it  is  less  than  the  middle  slope. 
The  program  has  the  ability  to  adjust  the  arbitrary  values  of  time  to 
ensure  the  middle  slope  is  the  largest.   Once  this  basic  shape  is  obtained 
the  program  calculates  the  value  of  time  at  which  the  parabola  has  a  maximum. 
This  process  is  repeated  using  the  calculated  maximum  time  and  slope  as  the 
middle  point  and  the  two  nearest  points  from  the  preceeding  iteration. 
When  the  difference  of  the  new  value  of  maximum  slope  and  the  previous 
calculated  maximum  slope  is  less  than  EPl,  a  program  input,  the  slope  is 
accepted  as  the  maximum.   Then  the  values  of  time,  slope  and  input  para- 
meters are  printed  out. 
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TABLR  2.      RESULTS  OF  SLO3  AND  SL05  FOR   7\    ■  0,0 
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TABLE  ^.     RESULTS  OF  SL02  AND  SLQ*f  FUR  g  =  0.0 
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6.    Discussion  of  Results. 

Table  1  indicates  that  programs  SL03  and  SL05  are  not  sensitive  to 
G  and  that  both  programs  give  the  same  solution  to  six  decimal  places 
for  NTU  greater  than  5.   The  variation  of  G  from  -0.1  to  +0.1  is  much 
larger  than  the  variations  used  in  the  SL02  and  SL04  programs  which  are 
more  sensitive  to  G.   All  results  presented  were  calculated  with  eight 
point  formulas. 

One  way  to  increase  the  accuracy  of  the  program  is  to  increase  the 
number  of  terms  in  the  partial  sum.   There  is  a  test  in  the  program 

which  will  increment  the  number  of  partial  sums  if  the  nth  average  has 

-9 
not  converged  to  within  1  x  10   of  the  preceeding  nth  average  for  each 

calculation  of  slope.   The  test  for  the  sufficient  number  of  partial 
sums  indicates  that  for  NTU  greater  than  five,  twenty  partial  sums  will 
give  the  same  accuracy  as  100  partial  sums  for  programs  SL03  and  SL05. 
At  NTU  of  five,  the  SL03  program  automatically  increased  the  number  of 
partial  sums  to  25  which  then  gave  as  accurate  a  solution  as  95  partial 
sums.   There  is  an  indication  that  for  lower  values  of  NTU,  more  partial 
sums  may  be  necessary. 

R.  E,  Maxim  [ioj  presented  a  table  of  values  of  maximum  slopes  cal- 
culated with  Schumann's  solution  for  the  special  case  of  ri   «  0,0. 
Schumann's  solution  is  an  infinite  series  of  Bessel  functions.   Each 
value  of  slope  was  calculated  using  Bessel  functions  accurate  to  ten 
significant  places  and  the  summation  was  continued  until  there  was  no 
change  in  the  tenth  decimal  place  of  the  solution.   Unfortunately s  the 
search  by  Maxim  could  not  distinguish  between  relative  maximum  slope 
values  and  his  search  ended  when  the  slope  at  the  next  increment  of  time 
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was  less  than  the  preceeding  slope-   Relative  maximum  slopes  have  been 
found  by  the  theoretical  solution  presented  in  this  paper  and  their  ex- 
istence is  supported  by  experimental  evidence  found  by  J.  M.  Bannon  ill  , 
Figure  26.   The  results  of  Table  2  clearly  indicate  that  Maxim's  values, 
though  accurate  to  six  places,  were  calculated  at  a  time  much  too  low  to 
have  found  the  largest  relative  maximum  slope  in  the  NTU  range  less  than 
five.   For  values  of  NTU  equal  to  and  greater  than  five,  both  the  time 
at  which  maximum  slope  occurred  and  the  values  of  maximum  slope  times  NTU 
calculated  by  SL03  agree  to  at  least  seven  decimal  places  with  Maxim's 
results.   SL05  solutions  were  in  complete  agreement  with  the  results  of 
SL03. 

The  computation  time  depends  on  the  degree  of  accuracy  required  to 

calculate  one  value  of  slope.   The  accuracy  required  of  the  search  routine 

-9 
for  the  results  of  Table  2  was  +  1  x  10   difference  in  the  new  value  of 

maximum  slope  as  compared  to  the  previous  iteration.   For  values  of  NTU 
less  than  ten,  the  number  of  iterations  required  was  four  or  five  times 
the  number  for  NTU  values  greater  than  ten.   Once  the  approximate  value 
of  time  at  maximum  slope  is  known,  the  search  routine  can  be  set  to  find 
the  maximum  with  fewer  iterations.   But,  even  without  optimum  starting 
times  the  computation  time  was  approximately  half  a  minute  for  the  7\  =   0.0 
case.   For  non  zero  values  of   /\       ,  the  computation  time  per  search  is 
expected  to  be  twelve  minutes. 

Results  of  SL05  for    /^    =.04  are  shown  to  be  insensitive  to  G  in 
Table  3.   One  test  point  at  NTU  =  10  indicates  good  agreement  between 
SL05  and  SL03  programs.   As  for  the  pi    -   0  case,  results  indicate  good 
agreement  with  C.  P.  Howard's  results  [c]  .   For  NTU  greater  than  five 
Howard's  results  are  consistantly  low  but  all  have  less  than  1%  differ- 
ence as  compared  to  SL05.   As  in  the  f\    =  0  case,  Howard's  results  for 
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NTU  equal  to  two  are  much  lower  than  those  obtained  by  SL05.   Again  this 
is  probably  a  failure  of  his  search  technique  to  distinguish  between 
relative  maximum  slopes.   Unfortunately,  he  did  not  publish  the  times 
at  which  the  maximum  slopes  occurred.   Howard  also  used  the  incremental 
step  search  and  the  search  ended  when  the  newly  calculated  slope  was 
less  than  the  preceeding  slope. 

Programs  SL02  and  SL04  are  quite  sensitive  to  parameter  G.   The 
results  obtained  graphically  and  listed  in  Table  4  were  compared  to  the 
accurate  solution  of  Maxim  for  f\    =0.0.   The  percent  error  of  the  graphi- 
cal results  was  less  than  +  1.5%  except  at  NTU  =  2.   Figure  6  indicates 
that  this  intercept  was  not  determined.   Considerable  effort  would  be 
required  to  find  this  intercept  and  obviously  the  mean  value  at  NTU  =  2 
has  unacceptably  large  error.   For  the  other  values  of  NTU,  however, 
intercepts  were  easily  found  as  indicated  by  Figure  5. 

Table  5  shows  similar  results  for  /\     =  0.005.   No  accurately  known 
values  are  available  for  this  region.   C.  P.  Howard's  results  are  com- 
pared to  the  results  obtained  by  SL02  and  SL04.   If  as  expected,  the 
SL02  and  SL04  results  are  1  to  1.5%  too  large,  then  C.  P.  Howard's 
solutions  for  /\      =  0.005  would  be  approximately  .5%  too  low  for  NTU 
greater  than  10.   Figure  7  again  shows  how  the  intercept  or  best  G  was 
obtained. 
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7.   Conclusions  and  Recommendations. 

Program  SL05,  using  the  VR,  Gauss -Legendre  formulas,  is  a  fast 
accurate  tool  capable  of  supplying  all  the  theoretical  maximum  slope 
data  needed  for  experimental  testing  of  compact  heat  exchanger  surfaces. 
The  SL03,  VI,  program  is  equally  accurate  but  slightly  slower.   The  SL05 
program  is  being  applied  to  verify  the  NTU  versus  Maximum  Slope  curves 
presented  by  C.  P.  Howard  [5j  .   The  results  of  these  calculations  will 
be  published  at  a  later  date. 

The  combination  of  Laplace  transforms  with  numerical  inversion  should 
be  applicable  to  most  boundary  value  problems  in  which  one  is  able  to 
express  the  solution  as  an  inverse  Laplace  transform.   Therefore,  this 
technique  is  primarily  limited  in  application  to  linear  systems  of  part- 
ial differential  equations  with  constant  coefficients. 

The  two  main  problems  encountered  in  the  numerical  inversion  were 
the  slow  convergence  of  the  integration  and  radical  behavior  of  the  func- 
tion of  s  at  time  equal  zero.   The  VR  and  VI  functions  were  cyclic  in 
nature  with  slowly  decreasing  magnitudes  and  period.   For  this  particular 
solution  the  Legendre  polynomial  approximated  the  functions  better  than 
Chebyshev's  polynomial.   However,  for  other  functions  a  number  of  other 
polynomials  may  give  better  results.   The  acceleration  of  convergence  is 
adequately  accomplished  by  the  present  numerical  inversion  but  it  gener- 
ates a  problem  in  addition.   Because  of  the  averaging  technique  applied 
to  the  partial  sums  to  accelerate  convergence  the  formal  error  analysis 
normally  applied  to  each  partial  sum  is  impossible.   Therefore,  the  only 
way  remaining  to  determine  the  accuracy  of  the  calculated  maximum  slopes 
was  by  comparison  with  R.  E.  Maxim's  accurately  determined  values  for  the 
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limiting  case  of  /{    =  0.   Even  though  direct  inversion  of  the  Laplace 
transform  was  possible  for  this  case,  the  numerical  inversion  was  ac- 
complished in  seconds  with  the  same  accuracy  attained  by  the  approxi- 
mation to  an  infinite  sum  of  Bessel  functions  carried  out  to  six  decimal 
place  accuracy.   For  non-zero  s\    >  there  were  no  accurate  solutions  with 
which  to  compare.   For  these  values  an  indication  of  accuracy  was  obtain- 
ed by  calculating  the  same  points  with  an  independent  program.   The  SL03 
program  is  not  independent  of  SL05  in  a  strict  interpretation  of  indepen- 
dence.  However,  they  are  different  programs,  SL03  using  the  imaginary 
part  of  the  inversion  integral  and  SL05  using  the  real  part.   Either 
program  should  give  the  same  result.   The  same  result  to  six  or  seven 
decimal  places  was  obtained  with  these  programs  indicating  excellent 
accuracy  for  a  numerical  process.   The  comparison  of  SL05  results  with 
Howard's  results,  obtained  with  a  rather  crude  finite-difference  tech- 
nique, indicated  that  he  obtained  very  accurate  results  by  judicious  ap- 
plication of  the  finite-difference  technique.   Most  of  his  results  were 
less  than  one  half  of  one  percent  low  when  compared  to  SL05. 

The  following  recommendations  indicate  avenues  of  investigation 
that  could  be  undertaken  with  either  computer  programs  developed  in  con- 
junction with  this  thesis  or  the  mathematical  technique  presented. 

It  is  recommended  that  the  temperature  programs  TEMP4  and  TEMP2  be 
applied  to  obtain  a  time- temperature  history  for  comparison  with  corres- 
ponding experimental  results  to  determine  how  well  the  mathematical  model 
represents  the  behavior  of  the  experimental  test  rig. 

An  investigation  should  be  conducted  to  determine  the  cause  of  and 
physical  significance  of  the  relative  maximum  slopes  obtained  both  experi- 
mentally and  theoretically  for  the  heating  cycle.   One  objective  would  be 
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to  determine  if  the  relative  maximum  obtained  experimentally  correspond 
in  magnitude  and  time  of  occurance  with  those  obtained  by  theory,   Ex- 
perimentally it  was  noted  that  the  relative  maximums  are  practically 
indistinguishable  for  the  cooling  of  the  solid,  but  that  for  the  same 
matrix  and  mass  flow  of  fluid  the  relative  maximums  are  very  distinct 
for  the  case  of  heating  the  solid,   An  investigation  should  be  made  to 
determine  theoretically  the  effect  of  longitudinal  conduction  for  the 
cooling  cycle  as  compared  to  the  present  solution  for  the  heating  cycle. 
The  objective  of  this  investigation  would  be  to  determine  if  the  theoreti- 
cal solution  is  capable  of  predicting  a  difference  in  response  on  cool- 
ing as  compared  to  heating.   The  cooling  cycle  solution  would  require 
appropriate  sign  changes  in  the  heat  balance  from  which  governing  partial 
differential  equations  were  derived.   This  may,  therefore,  require  a  new 
solution  for  cooling  since  the  sign  changes  could  cause  a  change  in  the 
auxiliary  cubic  equation  3b  of  Section  3.   It  is  possible  that  a  sign 
change  of  a  term  in  this  equation  could  cause  a  change  in  the  significance 
of  the  parameter  p\    on  the  result  of  fluid  temperature  or  slope  of  fluid 
temperature.   Bannon  {_1J  has  noted  experimentally  that  the  largest  maximum 
slope  obtained  by  cooling  agreed  fairly  well  with  that  obtained  by  heating. 

The  theoretical  results,  in  the  region  of  NTU  =  2,  are  difficult  to 
obtain,  and  the  determination  of  NTU  given  an  experimental  maximum  slope 
in  this  NTU  region  is  impossible ,   An  attempt  should  be  made  to  develop 
a  curve  matching  technique  employing  a  least  square  error  method.   This 
method  would  use  three  or  more  values  of  time  and  their  corresponding 
temperatures  from  an  experimental  time-temperature  history.   The  TEMP4 
program  would  calculate  theoretical  values  of  temperature  at  the  given 
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times  for  an  estimated  value  of  NTU.   Then  an  NTU  search  routine  could 
calculate  the  sum  of  the  squared  error  of  each  temperature.   After  re- 
peating this  operation  for  two  bracketing  values  of  NTUS  a  Lagrange 
iteration  method  would  calculate  the  NTU  with  the  least  squared  error. 
Another  possible  means  of  curve  matching  would  be  to  match  the  dimension- 
less  time  at  which  the  maximum  slope  occurs  for  the  region  of  NTU  =  2, 
rather  than  matching  the  magnitudes  of  maximum  slope.   This  would  require 
an  adequate  starting  mark  on  the  experimental  time-temperature  history 
and  a  correction  of  the  time  at  the  maximum  slope  to  account  for  the 
system  time  delay  of  the  experimental  response. 

A  more  general  system  of  governing  partial  differential  equations 
can  be  obtained  from  the  heat  balance  as  derived  b-*   Creswick,  by  in- 
cluding a  term  to  represent  the  energy  stored  in  the  fluid.   This  term 
is  required  for  transient  heat  transfer  testing  using  a  liquid  rather 
than  a  gas.   The  result  of  this  addition  changes  equation  2  of  section 
3  to  the  following 


dw  -     i 


3X     MTU 


u-nr-^f3Ar 


9t 


J    ? 


where   /^  ^  £.  M^  Af  **-*        ,  the  dimensionless  fluid  capacitance, 

Ws  Ca  3 

and   "^   is  the  density  of  the  fluid  in  lb/ft.   This  equation  combined 

with  Equation  1  of  Section  3  would  then  be  the  new  governing  partial 

differential  equations.   The  present  mathematical  technique  should  be 

capable  of  solving  this  system  of  equations. 
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APPENDIX  I 

Solution  for  Zero  and  Infinite  Longitudinal  Conduction. 
1.   ZERO  LONGITUDINAL  CONDUCTIVITY 

As  in  the  solution  for  finite  conductivity  the  development  of  this 
case  is  the  same  initally;  cf,  Equation  3a,  page  8. 

(i)    anr Cx,s)     h  aV(Xis)   b(s+i)jtor(x,s)    b25 nr(x^=n 
<3x3  a*2"  a.         ax  a 

For  this  case  where  /^  is  defined  as  a,  it  is  necessary  to  multiply  by  a 
The  substitution  of  a  =  0  results  in  the  equation 

(2)     or'  +  b  S-  nr^o 

S+l 

The  general  solution  of  equation  2  is 


(3)      nr(x,s)  «    C.(s)eh,y=  Ci««e"m< 


The  boundary  conditions  are  the  same  as  in  the  previous  solution. 
Boundary  conditions  3.  and  b  were  applied  in  order  to  arrive  at  equation 
1,  above.   Boundary  condition  C  is  v(0,t)  =  C  which  transforms  into  v(0,s) 
C/s.   Applying  boundary  condition  0  to  equation  3  results  in 
AT(0,  S)=  —  =  C|(S)  ?  therefore 


(4)    nr(x,s)=  c  0-  kl*      and  /T/Yi,s;  =  C_e    s^r 


It  follows  directly, 
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(5)  at       ~  Le 


and  ^t 


=  fe-£ 


This  special  case  was  solved  analytically  in  1929  by  T.  E.  W. 
Schumann  jl4j  •   The  following  development  will  show  that  the  present 
solution  agrees  with  Schumann's. 


Then 


Let  s'  =  s  +  1  and  apply  it  to  equation  4. 

_b  J~i  X 
AT(X,  5-1)  -_£__  C     5' 


or 


<6>    Ar(x,S-l)  =,  Ce"^^ 


w-i 


From  page  229  of  reference  [2 J  ,  we  find 


(7)  p-*r±  * 


X 


)         W 


/* 


.j.Cey/KT). 


Comparing  equation  6  and  7,  the  inverse  transform  of  equation  6  becomes 


From  page  294  Qf]  is 

£l{Hi-0-)}  =   e**-^"1^^)]  therefore, 


CD         -_      W 


(9a)  nrcxj)  --Ce~h(>"VE  (£)  *■  lM(zJZtt) . 


sn-zo 
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Define  a  =  btx  and  W  =  2  -ya   ,  then  compare  equation  9a  with  the  form 


w>  a"  2  ImUVa?). 


Note  that 


/t\r        t  ^^  -2 


Equation  9a  then  becomes 


do,  ortM-Ce         £ tma?  l^(zV^). 


Schumann's  solution  is 


(id   arts,*)  =  e"*-*^  z"  dlJ^(^^-V^) 

Equations  10  and  11  are  equivalent  which  can  be  shown  as  follows.   From 
page  392  [_4J  t    the  modified  Bessel  function  identify  is 

dWL  J  where   m±  O 

Then,    for  n  =  0  and  W  =  2  -]/&] 

dw  A 

Apply  the  chain  rule  such  that 

dCY]    =    dY      ^V 

d  a  dvv  '  da- 

it   follows  that 


djx0  c^Tio?)]  .  a  al,faVo:), 

da 
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By  mathematical  induction  it  can  be  shown  that 

(12)   d^.  [l»  (2-^)J  =  a^lm(2{^). 

dcT 

Equation  10  then  becomes 

as)  /vwt)«Ce-bx+tft*H!Cr.f2Va3 

From  page  392  [4]  , 

i»^i  =  i'" xjn(zi -f5?)y  and 

Note  that  Schumann's  development  is  for  the  case  of  heating  the  fluid  as 
time  increases  and  that  the  present  case  is  developed  for  cooling  of  the 
fluid  with  increasing  time,  which  changes  the  sign  of  t  in  the  exponential 
term         It  then  follows  that  equation  13  is,  indeed, 

Since  the  program  for  numerical  inversion  of  equations  5  and  6  is  already 
available  it  is  much  easier  and  probably  much  faster  to  apply  it  than  to 
evaluate  equation  11.   The  same  theoretical  solution  was  obtained  by  G.  L. 
Locke  in  1950  [9 J  . 


43 
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2.   INFINITE  LONGITUDINAL  CONDUCTIVITY 

When  f\     =  00  >  the  solid  temperature  is  no  longer  a  function  of  x. 

^  ^  N  (     ±\ 
In  fact  the  following  is  true:  -^ —  ^  >W-  0.     Therefore,  equation 

1  ,  page  7,  becomes  indeterminate  and  must  be  replaced  by: 

(1)  <*}±tt)_     C      rf\r(xft)-lX(i)\   dX  and  equation  2  is 
dt                  Jo      l                              -1 

(2)  2nHKt)    -  NTu[u.(t)-AT(X)t)] 

The  boundary  conditions  to  be  applied  to  this  case  are: 

a.  v(x,0)  =  0, 

b.  v(x,0)  =  0,  and 

c.  v(0,s)  =  C/s. 

Solve  equation  2  for  u(t)  and  take  the  derivative  of  u  with  respect 
to  t.   Then  substituting  these  back  into  equation  1  and  applying  boundary 
condition  c,  results  in 

MTU    2x2t  dt  ^    L  NTU3X  -lu* 

=  =±-  IV(l,t)-Cl  . 

NTU  L  -J- 

Multiply  through  by  NTU  and  rearrange  the  terms  to  arrive  at 


(3)  l^x^K  MTu^r(x^) ,  ar(i,t;-C=  o. 


Transformation  of  each  term  in  3  by  Laplace  gives 
For  further  discussion,  see  page  47, 
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(4)   anrlx.o)  >S<2nr(x,S)   .  MTuf-nr^x;o)  +  Snr(x,5)~]  +nr(i,s)-Q  -o 
^      PX        L  J        s 


ay 


Applying  boundary  conditions  a  and  b  and  dividing  by  s,  equation  4  reduces 
to 


3*  S^  ~5    # 


Again  by  treating  s  as  a  pa rameter^ equation  5  is  a  total  differential 
equation.   The  solution  of  the  differential  equation  is 


W   AT(X,S),  CxCs^e"1-   + 


-(ntu)x 


NTU   '  S 


Now  apply  boundary  condition  c  to  equation  6  which  results  in 


(7) 


'    S  NTU   S 


"2. -was) 


3 


Therefore,       Ci 


»•%[*- 


1 


NTU  •  SJ 


+ 


•1  _      /v(i,S) 
NTU  *S 


and 


(8) 


-NTU  «X 


1 


/V(x;S)  =  ^.e 

;  S  NTU   «S 


C    .AT  (IS) 
L  S 


-MTU"«X\ 


Solve  equation  8    for  v(l,s     as    follows: 


nr(l,S)=ke  + 


NTU  -sLS 


or(X^) 


_i 


Ci-e 


siTU      \ 


nr(l^) 


1  + 


(l-e 


^TU     j 


NTU    •   S 


0 

s 


(±-e'NTU  J  ^e-NTU 


NTU     .  $ 
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/  -^ru\  -MTU 

Define  £        _        (  1  -  €>         J       ;      then      £  -  1  ~  V  NTU  ,    and 

J      ~"  Niru 

The  particular  fluid  temperature  of  interest  is  that  temperature  at  x  =  1. 
Also  the  value  of  C  must  be  dimensionless,  for  comparison  with  Hondt's 
solution,  and  can  be  set  equal  to  one.   It  follows  directly  that 


fW(XS)    =  1  (1+  |  -  f'NTll) 

(s+j) 


(10)     nrd,s)  =       _i 1_niu_  / 

5+/       S^  S(S+?) 


The  inverse  Laplace  transform  of  equation  10  is 


-ft  _  ft  ftt 


which  reduces  to 
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(id    rril.t)  =  ±-  mtu  ye        , 


a      /-e"Nru 

where  \    -  . ,   The  solution  as  a  function  of 

N  T  U 


time  is  p-"™)  , 

This  solution  is  in  complete  agreement  with  that  derived  by  J.  R, 
Mondt  [if]  in  1961 .   Since  a  convenient  inverse  transform  was  available 
it  was  unnecessary  to  apply  the  intricate  numerical  inversion. 

To  derive  Equation  (1)  for  this  special  case  it  was  necessary  to 
refer  to  the  heat  balance  from  which  Creswick  j~3j  derived  the  governing 
equations.   The  last  term  of  Creswick' s  Equation  (1)  was  eliminated  and 
the  remaining  terms  were  arranged  in  terms  of  dimensionless  parameters. 
Then  both  sides  of  the  equation  were  integrated  over  the  length  of  the 
matrix  to  arrive  at  Equation  (1)  on  page  44. 
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APPENDIX  II 
Flow  diagrams  and  program  listings. 

The  flow  diagrams  will  preceed  the  applicable  program  listings 
The  symbols  used  in  the  flow  diagrams  are  defined  as  follows: 


(  ) 


o 


a 


Terminal   -   the  beginning,  end,  or  point  of 
interuption  in  the  program 


Input/Output   -   input  symbols  are  preceeded  by 
a  terminal. 


Decision  -  branches  are  labeled  by  sign 
according  to  the  sign  of  the 
value  enclosed. 

Connector  -  enclosed  numbers  refer  to  the 
prefix  numbers  of  the  corres- 
ponding statements  on  the  pro- 
gram listing. 

Offpage  Connector   -  designates  an  entry  to  or  exit 

from  a  page. 


There  are  six  distinct  programs  which  are  used  separately,  although 
they  could  be  combined  into  a  larger  package  to  make  use  of  one  or  another 
at  the  user's  option.   There  are  the  four  "slope"  programs,  SL02,  SL03,  SL04 
and  SL05  which  have  been  described  in  detail  in  the  body  of  the  thesis,  and 
the  two  "temperature"  programs  TEMP2  and  TEMP4,  which  have  been  referred  to 
in  the  text  but  not  described  in  detail.   The  latter  are  intended  to  pro- 
vide- a  time -temperature  history;  no  use  of  these  programs  has  been  made  in 
this  thesis.   See  page  35  for  recommendations  concerning  their  use. 


48 


PROGRAM  SL03~ 


(     Input       y~ Qh 


G,   B,   A,   It, 

TM,   EP1 ,    EP, 
L1,   LD,    11 


Set  11=0 
on  last 
data  card 


(~  Sto-p   10 


(Output) 
Program  isj 
Slope  3 


(Output) 
Evaluation/ 
Subroutine/ 
\2   Short 


(Output) 
Uses  Legendre 
Polynomial 


-<y 


Do  203  J=1,100. 
Call  Gauss Q 


*Also  applies  to  program  SL05 


€> 


L=L+LD 
Call  Aver" 


L 


Continue 


202 


■<* 


(Output) 
T1,   Slope 


(Output) 
G,   B,A,TT 

i SUM,    DIF, 
\     SUCJ) 


Call   SLOMAX 
!TM1=8TM/3 


(Output), 


't^: 
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PROGRAM  SL03 

DIMENSION  VR(2C),VI(20) ,SU(200),AV( 1  2  ,  1  0  ) ,  R  (  1  29  )  ,  X  (  1  29) ,ER(2C0) 
1RR(128),RI(123),T(10),SLC(  10) 
COMMON  RR.RI ,T,SLOfGt Rt A,PI,SLOPEfMT,NTl tNT2tEPl , II ,T1 
COMMON  NT3,TM 

101  REAC  100,G,B,A,T1.TM,EP1,EP,L1,LD,I 1 
100  F0RMAT(7F10.5,3I3) 

C      SET  11=0   ON  LAST  DATA  CARD         ,' 
C      II  IS  THE  AVERAGE  AND  THE  POLYNOMIAL  ORDER 
IF( I  1 )10, 10, 102 

102  CONTINUE 
NT  =  -1 

.  NT2=0 
K=l 

NT1=-1 

PI=3.  1U15926536 
WRITE  OUTPUT  TAPE  4,301 

301  F0RMAT(2UH  PROGRAM  IS  SLOPE  1,'N=8//) 
WRITE  OUTPUT  TAPE  U.302 

302  FORMAT(30H  EVALUATION  SUBROUTINE  2  SHORT/) 
WRITE  OUTPUT  TAPE  U,305 

305  F0RMAT(25H  USES  LEGENDRE  POLYNOMIAL) 
204  L=L1 

NT3=-1 

SU  =  0. 

D  =  0.0 

E=1  .0 

DO  203  J=1 ,100 

CALL  GAUSSC(D,E,VAL) 

D=U+1.0 

E=E+1.0 

SU( J)=VAL+SU( J-l ) 
503  IF( J-L)203,206,206 
206  L=L+LD 

CALL  AVERt SU,AV, 1 1  ,  J  ,  DI F , SUM) 
500  IF(DIF-EP)202,202,203 
203  CONTINUE 

CALL  AVER(SU,AV,  8 , 1 00, CI F , SUM ) 
202  EG=EXPF(G*T1 ) 

EG2=2.*EG/T1 

SL0PE=EG2*SUM 

PRINT  3,G,B,A,T1 

PRINT  i»,SUM,DIF,J,SU(  J) 

PRINT  2, T1, SLOPE 

221  CALL  SLOMAX 
TM1=8.*TM/3. 
IF(NT2)M06.U06.101 

U06  IF(T1-TM1)20M,20U,222 

222  PRINT  3,G,R,A.T1 

2  F0RMAT(3H  T=F10.5,3X,6HSL0PE=E20.8M 

3  F0RMAT13H  G= Fl 0 .5, 3X ,2HE=F 10. 5, 3X,2HA=F 10. 5 , 3X, 2HT=F 1 0. 5/ ) 

4  FCRMAT(5H  SUM=E20.8,3X,UH0IF=F10.9, 3X,3HSU( I3,2H)=E20.3/) 
GO  TO  101 

10  STOP  10 
END 
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PROGRAM  SL05 

DIMENSION  VR(2C),VI(20) ,SU(200),AV( 1  2  , 1  0  ) ,  R  {  1  29)  ,  X(  1  29) ,ER(200) 
1RR(128),RI( 128),T( 10),SLO(10) 
COMMON  RR.RI ,T t SLO ,G, B, A, P I , SLOPE ,N T, NT  1 , NT2, Epl , 1 1 ,T 1 
COMMON  NT3.TM 

101  READ  100,G,B,A,Tlf  TM , EP1 , EP, L 1 , LD,I  1 

100  FORMAT(7F10.5,3I3)  j 

SET  II  =  0   CN  LAST  DATA  CARD 
II  IS  THE  AVERAGE  AND  THE  POLYNOMIAL  ORDER 
IF(  11)10,10,102 

102  CONTINUE 
NT  =  -1 
NT2  =  0 
K=l 
NT1=-1 

PI  =  3. 1415926536 

WRITE  OUTPUT  TAPE  4,301 

301  F0RMAT(2UH  PROGRAM  IS  SLOPE  5,  N=8/ ) 
WRITE  OUTPUT  TAPE  4,302 

302  F0RMAT(30H  EVALUATION  SUBROUTINE  2  SHORT/) 
WRITE  OUTPUT  TAPE  4,305 

305  F0RMAT(25H  USES  LEGENDRE  POLYNOMIAL/) 
204  L=L1 

NT3=-1 

SU  =  0. 

D  =  0.0 

E=1.0 

DO  203  J=l ,100 

CALL  GAUSSX)(D,E,VAL) 

D=D+1.0 

E  =  E-H.0 

SU(J)=VAL+SU(J-1) 
503  IF( J-L)203,206,206 
206  L=L*LD 

CALL  AVER(SU,AV,I1.J,DIF,SUM) 
500  IF(DIF-EP)202,202,203 
203  CONTINUE 

CALL  AVER1SU.AV,  8, 1 00, DIF ,SUM ) 
202  EG=EXPF(G*T1 ) 

EG2=2.*EG/T1 

SLCPE=EG2«SUM 

PRINT  3,G,B, A,T1 

PRINT  U,SUM,DIF,J,SU( J) 

PRINT  2,T1 , SLOPE 

221  CALL  SLOMAX 
TMl=8.*TM/3. 
IF(NT2)406.406.101 

406  IF(T1-TM1)204,204,222 

222  PRINT  3,G,B,A.T1 

2  F0RMAT(3H  T=F 10.5, 3X, 6HSL0PE=E20.8/ ) 

3  F0RMAT(3H  G=F10 .5, 3X,2HE  =  F 10. 5, 3X,2HA*F 10. 5 , 3X,2HT=F 10.  5/) 

4  F0RMAT(5H  SUM=E20.8, 3X, 4HD  IF=F1 0.9,  3X , 3HSU(  13, 2H ) =E20.8 / > 
GO  TO  101 

10  STOP  10 
END 
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PROGRAM  SL02 


Q     Input  ~~\: M 


or 


j 


G,  B,  A,  T1,  TM 
TD,  EP,  L1?  LD, 


11 


Set  11  =  0 
On  last 
data  card 


(Output) 
Program  is 
Slope  2 


DO  500,  1=1 
200 

500  SU(I)= 
0.0 


*Also  applies  to  SL04 


(Output) 
Uses 
Che by s he v 
Polynomiaj 


(Output) 
Evaluation 
Subroutine 
25 


-O- 


603  DO  200  1=1,  8 

T  =  PI*P 

CALL  S00TS2- 

CALL  EVAL2S 

VR(I)  =  VR 

VI(I)  =  VI 
200  P1  =  P1+C1./9.) 


€) 


Call  SLOMAX 
TM1  =  8TM 

3 


X  Output) 
\g,B,A, 
\  Tl 


(Output) 
SUM,  DIF 
202)-  -\  SU(J) 


(Output) 
T1,  Slope 


<& 
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PROGRAM  SLG2 

DIMENSION  RR(128).RI(  128) ,VI(20),V*( 20) ,SU(200),AV( 1 2,1  0 ) , R ( 129 ) , 
1X( 129) ,SLO( 10), T( 10) 
COMMON  SLOPE, SLO,T, NT, NT1,NT2,K,EP1  ,G,B,A, PI,T1,TM 
COMMON  NT3 

101  READ  100,G,B,A,T1,TM,TD,EP,L1,LD, II 
100  F0RMAT(7F10.5,3I3) 

SET  II  =  0   ON  LAST  DATA  CARD 
IF(  II  ) 10,10, 102 

102  C 1=1.0 

210  WRITE  OUTPUT  TAPE  4,301 

301  F0RMAT(24H  PROGRAM  IS  SLOPE  2,  N=4) 
WRITE  OUTPUT  TAPE  4,304 

304  F0RMAT(26H  USES  CHEBYSHEV  POLYNOMIAL) 
401  WRITE  OUTPUT  TAPE  4,302 

302  F0RMATI30H  EVALUATION  SUBROUTINE  2S      /) 
W1C=. 10776070/. 98480775 

W2C=. 08 3333 333/. 86602540 

W3C=. 045908434/. 64278761  . 

W4C=. 01 29975 31/. 342020 14 

NT  =  -1 

NT1=-1 

NT2  =  0 

K=l 

EP1=. 000001 

NT3=-1 
204  P=3.1415926536/T1 

DO  500  1=1,200 
500  SU(I)=0.0 

L  =  L1 

Pl=l./9. 

SGN=-1. 

SU=0.  • 

00  203  J=1 ,100 
603  00  200  1=1,8 

Y=P1*P 

CALL  S00TS2  ( G, B, A, Y , R, X, RR  ,  R I ) 

CALL  EVAL2S(RR,RI,G,B,Y,C1,VR,VI) 

VR(I)=VR 

VI(I)=VI 
200  Pl=PH(1./9.  ) 

605  SU( J)=SGN*(W1C*(VI (4)+VI(5) )+W2C* W  I ( 3) +VI ( 6 ) )+W3C*{VI(  2)fVI(7) )♦ 
1W4C*(VI< 1 )+VI(8) ) )*SU( J-l ) 

IF(J-l)  212,212,213 

212  PRINT  3,G,B,A,T1 

213  CONTINUE 
SGN=-SGN 
Pl=PU(2./9.  ) 

IF( J-L)203,206,206 
206  L=L+LD 

CALL  AVER(SU,AV,I1 .J,DIF,SUM) 

IF(DIF-EP)202,202,203 
203  CONTINUE 

CALL  AVER(SU. AV, 10, 1 00, DIF, SUM) 
202  EG=EXPF(G»T1 ) 

EG2=2.*EG/T1 

SL0PE=EG2»SUM 

PRINT  4,SUM,DIF,J,SU( J) 

PRINT  2, T1, SLOPE 

221  CALL  SLOMAX 
TM1=8.*TM/3. 
IF(NT2)406,406. 101 

406  IF(T1-TM1)204,204,222 

222  PRINT  3,G,B,A.T1 

2  F0RMAT(3H  T  =  Fl 0.5t 3X , 6HSL0PE  =  E1 5.8/ / ) 

3  F0RMATI3H  G  =  F 10. 5t 3X, 2HB=F 10.5, 3X .2 HA=F 10.5 ,3X.2HT=F 10.  5/ ) 

4  F0RMATJ5H  SUM*E1 5. 8, 3X, 4HDI F*Fl 0. 9,  3X, 3HSU (  I3f 2H ) *E 1 5.8 / ) 
GO  TO  101 

10  STOP  10 
END 


53 


PROGRAM  SL04 
v   DIMENSION  RRM28)  ,RI(  1  28  )  ,  V  I  (  20  ) ,  V*  (20)  ,SU(200),AV(  1 2 ,1  0) ,R( 1 29 ) , 
1X( 129),SL0( 10), T( 10) 
COMMON  SL0PE,SL0,T,NT,NT1, NT2,K,EP1  , G. B, A, PI , T  1  ,TM 
COMMON  NT3 

101  READ  100, G,B.A,T1,TM, TD.EP.Ll.LD.il 
100  F0RMAT(7F10.5,3I3) 

SET  11=0   ON  LAST  DATA  CARD 
IF( 11)10,10.102 

102  Cl  =  1.0  •>" 
210  WRITE  OUTPUT  TAPE  4,301 

301  F0RMAT(2UH  PROGRAM  IS  SLOPE  4,  N=4) 
WRITE  OUTPUT  TAPE  4,304 

304  F0RMAT(26H  USES  CHEBYSHEV  POLYNOMIAL) 
401  WRITE  OUTPUT  TAPE  4,302 

302  FORMAT(30H  EVALUATION  SUBROUTINE  2S      /) 
W1C=. 10776070/. 98480775 

W2C=. 083333333/. 86602540 

W3C=. 045908434/. 64278761 

W4C=. 01 2997531/. 342020 14 

NT  =  -1 

NT1=-1 

NT2=0 

K=l 

EP1=. 000001 

NT3=-1 
204  P=3.1415926536/T1 

DO  500  1=1,200 
500  SU( I)=0.0 

L=L1 

Pl=l./9. 

SGN=1. 

SU  =  0. 

DO  203  J=l, 100 

IF( J-l)600,600,601 

600  Pl=-7./18. 
GO  TO  603 

601  IF(J-2)602,602,603 

602  Pl=ll./18. 

603  DO  200  1=1,8 
Y=P1*P 

CALL  S00TS2  ( G, B, A.Y, R, X.RR.R  I  ) 
CALL  EVAL2S(RR,RI,G,B,Y,C1,VR,VI) 
VR( I )=VR 
VI (I  )=VI 
200  P1=Pl+( 1./9.  ) 

IF ( J-l )604,604,605 

604  SU(1  )=.5*SGN*(W1C*(VR(4)+VR(5))+W2C«(VR(3)+VR(6)  )+W3C*(  VR (2 ) +VR ( 7 ) 
1  )+W4C*(VR( 1 )+VR(8) ))+SU(J-l ) 

GO  TO  212 

605  SU( J)=SGN»(W1C»(VR(4)+VR(5) )+W2C» (VR ( 3) +VR ( 6 ) ) +W3C» ( VR(  2)+VR(7) )  + 
1W4C»(VR( 1)+VR(8) ) )+SU( J-1 ) 

IF(J-I)  212,212.213 

212  PRINT  3.G.B.A.T1 

213  CONTINUE 
SGN=-SGN 

Pl  =  Pl  +  (2./9.  ) 

IF( J-L)203,206,206 
206  L=L+LD 

CALL  AVERtSU.AV.Il.J.DIF.SUM) 

IF(DIF-EP)202,202,203 
203  CONTINUE 

CALL  AVER(SU,AV,10,100,DIF,SUM) 
202  EG=EXPF(G«T1  ) 

EG2=2.*EG/T1 

SL0PE=EG2»SUM 

PRINT  U,SUM,DIF,J,SU( J) 

PRINT  2, Tl, SLOPE 

221  CALL  SLOMAX 
TMl=8.«TM/3. 
IF(NT2)406,406,101 

406  IF(T1-TM1)204,204,222 

222  PRINT  3,G,B,A.T1 

2  F0RMAT(3H  T=F 10. 5, 3X , 6HSL0PE=E1 5. 8/ / ) 

3  F0RMAT(3H  G=F10.5, 3X, 2H8=F 10. 5. 3X,2HA=F 10. 5, 3X. 2HT=F 1 0. 5/ ) 

4  F0RMAT(5H  SUM=E15 . 8, 3X, 4HD IF=F10.9,  3X, 3HSU( 13, 2H ) =E 1 5.9 /) 
GO  TO  101 

10  STOP  10 
END 
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PROGRAM  TEMP1+* 


C     Input 


G,  B  A,  T,i 
TM,  TD,  EE, 
L1,  LD,  11 
1M1 


r    Stop  10  ") MO 


(Output) 
Program  is 
\ Tempera  - 
ture  h 


Do  203  J=1?  100 

Call  Gauss Q 


Continue 


(Output) 
.Evaluation 
\Subroutine 
1  Short 


>06 


-O- 


L=L+LD 
Call 
AVER 


r 


Set  11=0 
On  last  ' 
data  card 


{Output) 
Uses 

iLegendre 
'olynomial 


(Output) 
G.   B,,A,   T 


">vAlso  applies   to  program 
TEMP2 


®t 


T=T+TD 


£-TM 


-<3- 


(Output) 
T,   TEMP, 
SUM,    DIF, 

isu(J) 
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PROGRAM  TEMPU* 

DIMENSION  RR(128),RI( 12e ) ♦ V K 20 ) , V* ( 20 ) , SU( 2Q0) , AV( 1  2,1  0  )  ,R  (  1  29  )  , 
1X(129),ER(200) 
COMMON  G,B,A,T,TM,TD,EP,L1,LD,I1,M1  , RR, RI,PI 

101  READ  100,G,BtA,T,TM.TD,EP,Ll,L0,I 1,M1 
100  F0RMAT(7F10.5,3I3,I1) 

SET  11=0   CN  LAST  DATA  CARD 
IF (  1 1  )  10, 10, 102 

102  CONTINUE 

PI=3. 1415926536 

WRITE    OUTPUT    TAPE    4,300 
300    F0RMAT(30H    PROGRAM    IS    TEMPERATURE   '4,    N=8/) 

WRITE    OUTPUT    TAPE    4,302 
302    F0RMAT(30H    EVALUATION    SUBROUTINE    1    SHORT/) 

WRITE    OUTPUT    TAPE    4,305 
305    F0RMAT(25H    USES    LEGENDRE    POLYNOMIAL/) 
204    L  =  Ll 

L  =  L1 

SU  =  0. 

0  =  0.0 

E=1.0 

DO    203    J=1,100 

CALL    GAUSSQ(C,E,VAL,ER) 

0*0+1.0 

E*E+1.0 

SU(J)=VAL+SU(J-1) 

IF( J-L)203,206,206 
206    L=L+LD 

CALL    AVER(SU.AV»I1.JVDIFVSUM) 

IF(DIF-EP)202,202,203 
203    CONTINUE 

CALL    AVER(SU,AV,    8, 100, C IF , SUM ) 
202    EG=EXPF(G*T) 

EG2=2.*EG/T 

TEMP=    EG2*SUM 

WRITE    OUTPUT    TAPE    4,1.G,B,A,T 

1  F0RMATC3H    G=F10 .5, 3X.2HB=F 10.5, 3X.2 HA=F 10. 5. 3X,2HT=F 10.  5/ ) 
WRITE    OUTPUT    TAPE    4,2, T,TEMP, SUM. DI F , J, SU( J ) 

2  FORMAT! 3H    T=F10.5, 3X, 5HTEMP=E20.8,3 X,4HSUM=E20.8,3X,4H0 I F=F 1 0.9 , 
13X,3HSU(I3,2H)=E20.8////) 

T=T+TD 

IF(T-TM)204,204, 101 
GO    TO    101 
10    STOP    10 
END 


*TEMP2    is   the   same  except    statement    300. 
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Slomax  subroutine 


(      Input    ")- 


(Output) 
Slope  is 
zero 


33>^- 


T1  =4-1(2) 


39" 


T1=B+T1 


K=K+1 
NT=0 


C 


siope,  acr, 

J  NT1 ,  NT2, 
"1NT3,  EP1, 
\TM,  B 


(Output) 
Slope  is 
negative 


T1=TM 
K=K+1 
NT=1 


Return 
to  221 


X 


T1  =  iT(D 
NT  =  -T 
NT3  =  0 


Return 
to  221 


\(  Output) 
OUT  23 


(Output)/ 
K .  •  =  ■  V 


c 


Return 
to  221 


Return 
to  221 


J> 


\{  Output)) 
OUT  19 


/^Return-  ~^\ 
V  to  221    J 
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Slomax  Subroutine  (cont'd) 


SLOT 

=SL0lf 

|t(i)=iM 


(Output) 
OUT  25 


SL0^+=Slope 
T7=T(lO-T.(2) 


•DIFF  = 
SLOlf 
-SL05 


SL03=SL02 
T(3)=T(2) 
SL02=SL0lf 
T(2)=T(*0 
SLO?=SLQl4. 


(Output/ 
OUT   19 


NT2  =  1 


DIFF  = 
SLO*f 
-SL05 


Return        >j 
to   221        J 


SL01=SL02 
T(1)=T(2) 
SL02=SL0if 
T(2)=T(lf) 
SLOfeSLOk 


(Output; 

OUT   1 


SL03 

=SLOl+ 
T(3)=TCU) 


(Output) 
OUT  27 


1 
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41 

44 

42 

45 

43 

5 

2 


38 

39 
40 


19 


35 
21 


23 


36 
22 


30 


SUBR 

DINE 

CO.MM 

COMM 

IF(S 

WRIT 

FORM. 

GO  T 

WRIT 

FORN 

IF(N 

IF(N 

K=l 

T(K) 

SL01 

SLO( 

IF(N 

T  1  =  1 

GO  T 

T1  =  ( 

K=K  + 

NT  =  0 

RETU 

T(K) 

SL02 

SLC( 

SL05 

IF(S 

Tl  =  l 

NT=- 

NT3  = 

PR1N 

FORM 

RETU 

T1  =  T 

K=K  + 

NT=1 

RETU 

T(K) 

SL03 

IF(S 

Tl=4 

NT=1 

PRIN 

FORM 

RETU 

SLO( 

K=K  + 

C1=S 

C2  =  - 

C3  =  S 

C4=S 

C5=- 

C6  =  S 

Tl  =  . 

NT1  = 

PRIN 

FORM 


OUTINE  SLOMAX 

NSION  LTITLE(10),RR(12E),RI{123),T{10)fSLO(10) 

CN  RR.RI ,T,SL0,G,R,A,PI,SL0PE,^T,NT1,NT2,EP1  ,11 , Tl 

ON  NT3.TM 

L0PE)4l ,42, U3 

E  OUTPUT  TAPE  4,4U 

AT(20H  SLOPE  2  IS  NEGATIVE/) 

0  43 

E  OUTPUT  TAPE  4,45 

AT(26H  CALCULATED  SLOPE  IS  ZERD./) 

TD5.6,  6 

T)2,3,4 

=  T1 

=SLOPE 

K)=SLOPE 

T3)39,38,39 

.*T(2)/2. 

0  40 

B  +  TD/2. 

1 


RN 

=  T1 

=SLOPE 

K)=SLOPE 

=  SL02 

L02-SL01  )19,21,21 

.»T( 1)/4. 

1 

0 

T  35 

AT(7H  OUT  19/) 

RN 

M 

1 

RN 

=  T1 

=SLOPE 

L02-SL03)23,22,22 

.»T(3)/3. 


T  36 

AT(7 

RN 

K)  =  S 

1 

L01* 

SL02 

L03* 

L01* 

SL02 

L03* 

5*(C 

0 

T  30 

AT(3 


H  CUT  23/) 

LOPE 

(T(2)«T(2)-T(3)«T(3)  ) 
*(T( 1)»T(1 )-T{3)»T(3)  ) 
(T( 1 )*T(1 )-T(2)«T(2)) 
(T(2)-T(3)  ) 
♦  (T( 1)-T(3)) 
(T(l  )-T(2)  ) 
1+C2+C3)/(C4+C5+C6) 

»K 

H  K=I3/) 
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25 


31 

20 
8 


27 


32 

26 
11 


12 


33 

14 


15 


34 


RE  TURN 

ER1=0.C 

SLCU=SL 

T(K)=T1 

SL0(K)= 

T7=T(U) 

IF(T7)7 

DIFF=SL 

IF(DIFF 

SLG1=SL 

SL0(1)= 

T( 1 )=T( 

PRINT  i 

F0RPAT( 

GO  TO  1 

IF(DIFF 

CONTINU 

NT2=1 

RETURN 

OIFF=SL 

IF(DIFF 

SLC3=SL 

SL0<3)= 

T(3)=T( 

PRINT  3 

FORMAT! 

GO  TO  1 

IF(DIFF 

CONTINU 

NT2=1 

RETURN 

CONTINU 

SL01=SL 

SLC(  1  )  = 

T( 1 )=T( 

SL02=SL 

SL0(2)= 

T(2)=T( 

SLC5=SL 

PRINT  3 

FORPAT( 

GO  TO  1 

CONTINU 

NT2  =  1 

RETURN 

CONTINU 

SL03=SL 

SL0(3)= 

T(3)=T( 

SL02=SL 

SL0(2)= 

T(2)=T( 

SLC5=SL 

PRINT  3 

FORMAT( 

GO  TO  1 

END 


OPE 

SLOPE 

-T(2) 

,0,9 

0U-SLC5 

)25, 20,20 

OU 

SL04 

4) 

1 

7H  OUT  25/) 

-EP1  )  1U,1U,  15 
E 


0U-SLC5 

)27,26,26 

04 

SLOU 

U) 

2 

7H  OUT  27/) 

-EP1)  11,11,  12 
E 


E 

02 

SL02 

2) 

04 

SL04 

U) 

04 

3 

7H  OUT  12/) 


E 

02 

SL02 

2) 

04 

SLOU 

U) 

04 

4 

7H  OUT  15/) 
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33 


35 

J*  9 


37 


SUB 
DIM 

12  = 
N1  = 
DO 

A( 
Nl  = 

13  = 
Nl  = 
DO 
N1  = 
DO 
A(  J 
DIF 
SUN 
END 
SUB 
DIM 
COM 
COM 
IF( 
W=l 
XSI 
XS1 
XSI 
XSI 
XSI 
XSI 
XSI 
XSI 
XSI 
WSI 
WSI 
WSI 
WSI 
WSI 
WSI 
WSI 
WSI 
WSI 
Nl  = 

12  = 

13  = 
K2  = 
Al  = 
A2  = 
15  = 
00 
X(  I 
W(  I 
X(  I 
W{  I 
K2  = 
15  = 
SUM 
DO 
Xl  = 
CAL 
SUM 
RET 
END 


ROUT 
ENSI 
I  1  +  1 
Jl-I 

1  1  = 
III) 
Nl  +  1 
11-1 
12 

2  1  = 
N  1-1 
2  J= 
,1  +  1 
=  ABS 
=  A(N 

ROUT 

ENSI 

MCN 

MON 

W)5, 

.0 

(  16) 

(17) 

(18) 

(  19) 

(20) 

(21  ) 

(22) 

(23) 

(2U) 

(  16) 

(  17) 

(  18) 

(  19) 

(20) 

(21  ) 

(22) 

(23) 

(24) 

8 

19 

Nl/2 

12 

(E  +  D 

(E-D 

Nl 

35  I 

)=A1 

)=A2 

5)=A 

5)  =  W 

K2-1 

15-1 

=  0.0 

37  I 

X(I) 

L  EV 

=  SUM 

URN 


INE  AVER(S, A,  II. Jl,CIF,SUM) 
ON  S(200) ,A( 12, 1C) 

1-1 

1,12 

=  0.5«(S(N1 )+S(Nl+l )  ) 


1,13 

1,N1 

)=0.5»(A(J+1, I  )+A( J,I)  ) 

F( (A(N1 , 13+1) -A (Nl-1,  13+1)  )/A(Nl,  13+1)  ) 

1  ,13+1) 


,13+1) 

NE    GAUSSC(D,E,SUM) 

IN    XSI  (24),  WSI  (2i4),X(9),W(9)  ,RR(  128J.RK  128)  ,T(1  0)  ,SLO{  10) 
R,RI,TfSLO,GtBf AfPIt SLOPE,MT, NT1 ,NT2, EP1 , 1 1 ,T1 
IT3.TM 
,5 


=  0.  18 
=  0.52 
=  0.79 
=  0.96 
=  0.0 
=  0.32 
=  0.61 
=  0.83 
=  0.96 
=  0.36 
=  0.31 
=  0.22 
=  0.10 
=  0.33 
=  0.31 
=  0.26 
=  0.18 
=  0.08 


34346425 
55324099 
6666U774 
02098565 

4253M234 
3371U327 
6031 1073 
81602395 
2683783U 
37066459 
23810345 
12285363 
02393550 
23470770 
0610696U 
06481607 
127438814 


)»0.5 
)*0.5 

=  1,13 

-A2«XSI(K2) 

•WSKK2) 

1+A2»XSI(K2) 

(I) 


=  1,N1 

AL35(X1,VR,VI ) 
+VI*W(I) 


61 


SUBROUTINE  SCOTS2 ( R, X, Y ) 

DIMENSION  R( 129 ) , X ( 1 29 ) ,RR ( 128) .RI<  128) * ITC 10) 

COWCN  G,B,A,T,TM,  TO, EP,L 1 , LD, I  1 ,M1  ,RR,RI,PI 

R(l  )*1.0 

R(2)*B 

R(3)*-(B/A)«(6+1.) 

R{U)=-{ (B»«2)/A)»G 

X( 1 )=0.0 

X(2)*0.0 

X(3)=-(B/A)»Y 

XU)=-( (B*«2)/A)#Y 

N0  =  3 

MM  =  0 

IT  =  0 

CALL  TOOTS2(R.X,NO,IT,MM,-.5,*.5) 

IF(RR(1 ) )U,U,S 

RRR=0. 

RII=0. 

5  RRR  =  RR(  1  )  > 
RII-RK1) 

RR(1)=RR(2) 
RI(1)=RI(2) 
RR(2)=RR(3) 
RI (2)=RI(3) 
RR(3)=RRR 
RI(3)=RII 
GO  TO  3 
h    IF(RR(2))3,3,6 

6  RRR=RR(2) 
RII=RI(2) 
RR(2)=RR(3) 
RI  (2)=RI(3) 
RR(3)=RRR 
RI(3)=RII 

3  RETURN 
END 
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For  use  with  SL02  and  SL04. 


100 


102 
101 


SUBROUTINE  EVAL2S(RR,RI ,G,B,Y,C1,V*,VI) 

DIMENSION  Z(3),SR(4,U),SI{4,i4),E(3),S(3),C(3),RR(128),RI(128),W(3) 

DO  100  1=1,3 

Z(  I  )  =  ( (RR( I ) »RR( IJ-RI { I )»RI( I) )/B)  +  RRU ) 

W(  I  )=(2.»RR( I)»RI( I)/8)+RI(  I) 

D0101  1=1,2 

D0102  J=2,3 

SR(I,J)=(Z(I)*Z{J))-{W(I)«W(J)) 

SI ( I,J)=(Z( I )*W(J) )+(Z( J)#W( I ) ) 

CONTINUE 

E31=EXPF(RR(3)+RR( 1 ) ) 

S31=SINF(RI  (3)+RI(  1)  )  ' 

C31=C0SF(RI(3)+RI< 1) ) 

AR1=E31*C31 

AI1=E31«S31 

E3=EXPF(RR(3) ) 

S3=SINF(RI(3) ) 

C3=COSF(RI(3) ) 

DR1=E3*C3 

DI  1=E3*S3 

E1=EXPF(RR(  1  )) 

S1=SINF(RI( 1 ) ) 

CI l=COSF(RI ( 1) ) 

ER1=E1*C11 

EI  1  =  E1*S1 

SRN=SR(2,3)-SR( 1,2) 

SIN=SI(2,3)-SI( 1,2) 

SRD1=SR(2,3)-SR(  1,3) 

SIDl=SI(2,3)-SI(l,3) 

SRD2=SR( 1,3)-SR(1,2) 

SID2=SI( 1,5)-SI( 1,2) 

ANR=(SRN«AR1 )-(SIN»AIl ) 

ANI=(SIN*AR1 )+(SRM*AIl ) 

DNR1=SRD1*DR1-SID1»DI1 

DNI 1=SI01*DR1+SRD1«DI1 

0NR2=SRD2»ER1-SID2*EI1 

DNI2=SID2»lR1+SRD2*EI1 

DNR=DNR1+DNR2 

0NI=DNIUDNI2  > 

ANUR=ANR*DNR+ANI»ONI 

ANUI=ANI*DNR-ANR»DNI 

DENOM=DNR*DNR*DNI*DNI 

VR=C1»( ANUR/DENOM) 

VI=C1*(ANUI/DEN0M) 

END 
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100 


102 
101 


103 


Z?I)2((RrIi)»»2-RKI)»*2)/B)|RR(I) 
W(I)  =  (2.»RR(I)*Rnn/B)+RI(I) 
0101  1=1.2 

RII?J)  =  (Z(1I»Z(J))-IW(  I)«W(J)| 


D0101 
D0102 

s 
s 

CONTINUE 

E31=EXPF(RR(3)+RR  1) 

S31=SINF(RI(3)+RI 

C31=COSF(RI(3)+RI  1 

E32=EXPF(RR(3)+RR  2 

S32=SINF(RI(3)+RI  2 

C32=C0SF(RI<3|+RI  2 

E21=EXPF(RR(2)+RR    \\ > 

C21=COSF(RI(2)+RI     1 

S21=SINF(RI(2)+RK1>> 

DO    103    1  =  1 t 3  #  w  %  % 

E(  I  )=EXPF(RR(I ) ) 

S(  I)=SINF(RI(I) ) 

C(  I)=COSF(RI(I) ) 

AR1=E31*C31 

AI1=E31«S31 

AR2=E21»C21 

AI2=E21*S21 

BR1=AR2 

BI 1=AI2 
v BR2=E32»C32 
BI2=E32*S32 
CR1=6R2 
CI 1=BI2 
CR2=AR1 
CI2=AI1 
AR12=AR1-AR2 
AI 12  =  AI 1-AI2 
A1R=(SR(2,3) 
Al I=(SI(2,3) 
DR1=E(3)*C(3 
DI 1=E(3)»S(3 
DR2=E(2)«C(2 
DI2=E(2)«S(2 
DR12=DR1-DR2 
0I12=DI1-DI2 
ER1=E(1 )«C( 1 
EI1=E(1 )»S( 1 
ER2=DR1 
EI2=DI1 
ER12=ERl-ER2 
EI 12=EI 1-EI2 
FR1=DR2 
FI 1=DI2 
FR2=ER1 
FI2  =  E 1 1 
FR12=FR1-FR2 
FI 12=FI1-F12 
BR12=BR1-BR2 
BI12=BI 1-BI2 
CR12=CR1-CR2 
CI 12=CI1-CI2 
B2R=(SR( 1,3) 
B2I=(SI (1,3) 
C3R=(SR( 1,2) 
C  3  I  =  (  S  I  (  1 ,  2  ) 
D1R=(SR(2,3) 
Dl I = ( SI (2,3) 
E2R=(SR( 1," 
E2I=(S1 ( 1, 
F3R=(SR( 1, 
F3I=(SI ( 1,2) 
ANR=AlR+B2R+ 
ANI=A1I+B2I+ 
DNR*D1R+E2R+ 
DN1=D1H-E2I  + 
ANUR=ANR*ONR 
ANUR=ANR*ONR 
ANUI=ANI»DNR 
ANUI=ANI*DNR 
DENOM=DNR*DN 
DENOM=DNR»DN 
VR=C1*( ANUR/ 
VI=C1»(ANUI/ 
ENO 


For  use  with  SL02  and  SL04 
when  NTU  less  than  5 


•  AR12 

•  AR12 
) 

) 
) 
) 


)-(SI< 
)+(SR( 


2,3)*AI12) 
2,3)*AI12) 


3) 
3) 


*BR12 

•  BR  12 

•  CR12 
*CR12 

•  DR12 
*DR12 

•  ER12 
*ER12 

•  FR12 

•  FR12 
C3R 
C3I 
F3R 
F3I 

+  ANI« 
+  ANI* 
-ANR« 
-ANR* 
R  +  DNI 
R+DNI 
DENOM 
DENOM 


)-(SI( 
)+(SR( 
)-(SI( 
)+(SR( 
)-(SI( 
)+(SR( 
)-(SI( 
)+(SR( 
)-(SI( 
)+(SR( 


DNI 
DNI 
DNI 
DNI 

•  DNI 

•  DNI 
) 
) 


1,3 
1,3 
1,2 
1,2 
2,3 
2,3 
1,3 
1,3 
1,2 
1,2 


)«BI12) 
)*BI12) 
)»CI1?) 
)»CI12) 
)«DI12) 
)»DI12) 
)»EI 1 2) 
)»EI12) 
)»FI 1 2) 
MFI12) 
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For  use  with  TEMP2  anif  TEMP4 


SUBROUTINE    E VAL 53S ( X 1 , VR , V  I  ) 

DIMENSION    Z(3)  , SR(4,4),SI (  U  ,  U  )  ,  E  {  3 )  ,  S(3),C(3) ,RR( 128), RI (7  28),W(3) 
1 ,R( 129) tX( 129) 
COMMON    GfBtArTiTMt'TDf  EP'tLlf  LOt  MvNIf  RR»RItPI 
Y=(PI/T) *X1 
CALL    SOOTS2(R,X,Y) 
Cl=1.0 

DO    100    1=1,3 
Z(  I  )  =  ( (RR(  I )*RR(  I )-RI(  I  )*RI(  I  )  ) /B)  +  RR(  I  ) 

100  W(  I  )  =  {2.*RRM)*RI(  I  )/B)+RI  (  1  ) 
DO101     1=1,2* 

DO102    J=2,5 

SR(  I,J)=(Z(I  )*Z(J)  )-(W(  I)*MJ  )  ) 
102    SI ( I,J)=(Z( I )*W( J) )+( Z( J)*W( I)) 

101  CONTINUE 
F31=EXPF(RR(3)+RR( 1 ) ) 
S3",=SINF(RI(  3)+RI(  1  )  ) 
C31=C0SF(RI (3)+RI( 1 ) ) 
AR1=E31*C31. 
AI1=E31»S31 

E3=EXPF (RR(3) ) 
S3  =  SINF(RI ( 3)  ) 
C3  =  C0SF(RI (3)  ) 
DR1=E3*C3 
DI1=E3*S3 
E1=EXPF(RR(  1  )  ) 
S1=SINF(RI ( 1  )  ) 
CI l=COSF(RI(  1  )  ) 

ER1=E1*C11  •      • 

E  II  =  E 1  * S  1 

SRN=SR(2,3)-SR( 1,2) 
SIN=SI (2,3)-SI( 1,2) 
SRD1=SR(2,3)-SR( 1,3) 
S I D 1 =S I (2.3J-SI (1,3) 
SRD2=SR( 1 ,3)-SR( 1,2) 
SID2=SI( 1,3)-SI ( 1,2) 
ANR=(SRN*AR1)-(SIN*AI  1 ) 
ANI=(SIN*AR1 )+(SRN«AI7 J 
DNR1=SRD1*DR1-SID1*DI1 
'  DNI1=SID1 *DR1+SRD1  *DI1 
DNR2=SRD2*ER1-SID2*EI 1 
DNI2=SID2*ER1+SRD2*EI 1 
DNR=DNR1+DNR2 
DNI=DNI1+DNI2 
DENR=(G*DNR)-Y«ONI 
DENI=(Y*DNR)+G»DNI 
DENOM=(D£NR»OENR)+DENI*DENI 

VR=C1«C0SF<PI»X1 )*(ANR*CENR+ANI*DE\I )/DEN0M 
VI=-C1*SINF{PI*X1 )«( ANI*DENR-ANR*DENI  )/DENOM 
END 
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For  use  with  SL03*  and   SL05,      R    =  0. 

SUBROUTINE    EVAL35 ( X 1 . VR  ,  V  I  ) 

DINENSION    VR(20),VI(20)  ,SU(200) ,AV(  1  2  ,  1  0  ) ,  R  {  129  )  ,  X(  1  29)  ,ER(2CC) 
1RR{ 128) ,RI ( 128),T{ 10) ,  SL0( 10) 
C0KK0N    RR.RI  ,T,SLOfGf  B,  A, P I , SLOPE ,M T, NT  1 ,NT2 , EP1 , 1  1 , T 1 
COtfKCN    NT3,TM 
Y=(PI/T1)*X1 
Cl  =  l. 

01  =  ( (G*G)+G+(Y*Y) )/( <G+1. )«(G+1.)+(  Y*Y) ) 
E1=EXPF(-B«D1) 
F  =  Y/( (G+l. )*(G+1.)+(Y*Y)) 
C=CCSF(B*F) 
S=SINF(B*F) 

VR=C1»E1«C*C0SF(PI*X1  ) 
VI=C1*E1*S*SINF(PI«X1  ) 
END 


For  use  with  TEMP2  and  TEMP4,     A   =  0 

SUBROUTINE    EVAL3TC X 1 , VR, VI  ) 
^?RU29)I?X(129),SRU,U),  ^I(l4»U,'E(3,  tSC3),CC3),RR(128),RIC128),WC3) 
C0MM0NTd,B.AfT,TM,T0,EP,Ll,L0,Il,Ml  ,RR,RI,PI 

cT=i.o 

G1=G+1. 

DEN=G1*G1+Y»Y 

01=(G*G+Y«Y)/DEN 

F1=(G1»Y-G*Y)/DEN 

E1=EXPF(-B«D1) 

C=C0SF(B*F1) 

S=SINF(B*F1) 

DENCM=G*G+Y*Y 

yR=Cl»COSF(PI«Xl)#El«(C«G+S*Y)/DENaM 

VI=C1»SINF{PI«X1 )«E1«(S«G-C*Y)/DEN3M 

END 
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